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LIOUVILLE’S EQUATION AND THE n-BODY PROBLEM 


by 

R. G. Langebartel 
Goddard Space Flight Center 


SUMMARY 

The motion of a system of particles is examined on the basis 
of the fundamental equation in statistical mechanics. The Dirac 
delta function is used to describe systems which are discrete in 
position space, velocity space, or both as degenerate cases of 
continuous systems. The approximation procedure, necessitated 
by the nonlinearity of the problems, is based on the use of ex- 
pansions in successive derivatives of the delta function. This ap- 
proach leads to sum-difference- differential equations of a novel 
form for the quantities of interest, equations subject to a variety 
of techniques for solution. The method is applied to the dispens- 
ing and dispersion of the ’’West Ford Needles” belt, and to the 
problem of permanence of symmetry of the configuration of a 
collection of Echo-type satellites. 
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LIOUVILLE’S EQUATION AND THE n-BODY PROBLEM 


by 

R. G. Langebartel 
Goddard Space Flight Center 


INTRODUCTION 

The non-linearity of most problems in celestial mechanics usually necessitates the use of 
some approximation method. The scheme outlined in the present work leads to equations of a 
form rather different from those encountered in the standard procedures. The point of departure 
of the method is the regarding of the particle or collection of particles as constituting a medium 
in phase space subject to the fundamental collision-less equation of statistical mechanics, Liou- 
ville’s Equation. Such a medium is of necessity highly degenerate, so that its distribution function 
will involve the Dirac delta function. An example of a classical use of such ’’singularity functions" 
in statistical mechanics is the "microcanonical ensemble" of Gibbs for which the degeneracy is 
in the energy. However, for the celestial mechanics n-body problem the distribution function must 
be taken as degenerate in both velocity and position space. Some problems (e.g., the determination 
of the motion of the West Ford needles) call for a distribution function degenerate only in velocity 
space. A medium of this type is not the customary "gas" of statistical mechanics, but instead is 
the "completely incoherent medium" of Lichtenstein (Reference 1). 

The first sections are preliminary in character and introduce the use of the delta functions to 
describe the degenerate media. The approximation method is developed in the section devoted to 
the improvement of an approximate orbit. 

THE FUNDAMENTAL EQUATION 

The principal condition of the density distribution function is that it satisfy Liouville’s Equa- 
tion, the Equation of Continuity. This means that the density of the medium throughout its motion 
remains constant in phase space. The forces are assumed to be conservative, the force vector 
being the negative of the gradient of the potential. 

Let q r and p r be the n Hamiltonian coordinates and their conjugate momenta, and H(q r , p r , t) 
the Hamiltonian function. "Phase space" is the 2n-dimensional space of points 
(q x , • • • . q n , p lt • • * , p n ) . Let f (q k , p r , t) represent the density distribution function in phase 
space. Liouville’s Theorem (Reference 2, p. 83; Reference 3, p. 266) asserts that 

dj_ V /dH d f dH dj \ 

dt + /_\$P k dq k dq k dp k ) * (!) 

k 
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This first order partial differential equation in f is the fundamental equation for this theory. If 
the forces involved are purely external then the equation is linear. But if the self- gravitation of 
the particles is taken into account then f enters into h through the potential function v , and the 
equation becomes a nonlinear integro- differential equation which in almost every case must be 
handled by some approximation procedure. 

Liouville’s Equation (Equation 1) can be derived directly from the Hamiltonian Equations of 
Motion, but it is instructive for the case of a degenerate medium to show the equivalence of Equa- 
tion 1 to the Newtonian Theory by making use of the delta function S(x). Such a "function" has 
been employed for years but it is only relatively recently that a satisfactory theory incorporating 
them has evolved (References 4, 5, and 6). Physically S(x) represents the density distribution 
function of a unit mass concentrated at the origin. 

If the motion of a certain system with n degrees of freedom is 


= Qr K' Ps' 0 ■ 

Pr ^ P r K' Ps- t) > 

then the phase space density function for this system is 


(2) 


f = >K "Qr’ Pr~P r ) $ ‘''(qi-Ql. •••. Qn-On- Pl- P l’ •••• PH' P n) ■ (3) 

We want to show that substitution of Equation 3 into Equation 1 leads to a result equivalent to that 
predicted by the Newtonian Theory. With the notation 

di(q s -Q s , P s -P s ) 

S r = <5(q r -Q r ) 


t»(q 5 ~Q S . p s - p s ) 

5 n+ r ' (9(p r -P r ) 

Equation 1, for the substitution of Equation 3, takes the form 



( 4 ) 
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where H p ^ p = d 2 H (Q s , P s , t)/<9P r <9P k ; i.e., the original variables q r and p r in h are replaced by 
the functions Q r and p r . In deriving Equation 4 use is made of formulas for the multiplication of 
the delta function or its derivative by another function. The particular formulas employed in this 
instance are (Reference 5, p. 10) 


H qk K. p s , t)S r = \(Q s ,P s , t)8 r - (q s ,p 5 , t)S = H qt S r - H QkP S 


\ ( q s’ P s’ t ) S n+r H Q k 8 n+r H Q k P r 


5 . 


Setting the coefficients of S r and 8 n+r to zero gives a partial differential equation system in 
Q r and P r of a rather special form, a "quasi- linear system with same principal part" (Reference 7, 
V. II, p. 117). Its solution can be expressed as an arbitrary function of 2n functions of q r , p r , t. 
These functions in the present case turn out to be precisely the Newtonian integrals of motion for 
the given Hamiltonian. The coefficient of h in Equation 4 can be written in the form 



k 


and since g(x) b(x - a) = g(a) &(x - a) the b term is 



k 


^Pk 


b(q r -Q r , P r "P r ) 



k 


That Equation 4 represents the Newtonian Theory can be verified in a much briefer fashion 
by considering the special case Q r = Q r (t), P r - P r (t) , the form in which the solution to a prob- 
lem in dynamics is most often sought. For Q r and P r independent of q r and p r the vanishing of 
the coefficients of S r and 8 n+r gives immediately 


d Qr 

■dt = H P • 

r 

dP r 

"dF = " H g r : 


i.e., P r and Q r satisfy the Hamiltonian Equations of Motion. Note that for P r and Q r independent of 
p r and q r each term in the coefficient of S in Equation 4 vanishes. 

Other choices for the variables present in Q r and p r lead to certain specialized forms for the 
quasi-linear system derived from the coefficients of S r and S n+r in Equation 4. This permits a 
choice of equations to be used for a particular problem. 
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Example-Harmonic Oscillator 

The Hamiltonian is H = 1/2 (q 2 +p 2 ) and we will take Q = Q(t) and P = P(q) . From Equation 4 


__ dP 

p + q 


o , 


for which P = ± -|/c x 2 -q 2 . The remaining equation taken from Equation 4 is then merely 


dQ 

dt 




0 , 


and the complete solution is 


Q - c i sin (t - c 2 ) . 


P “ ± |/c J 2 - c J 2 sin 2 (t - c 2 ) - Cj cos (t - c 2 J . 

The equations taken from Equation 4 under the assumption P r = P r (t), Q r - Q r (p s , t) might 
be easier to handle than those resulting from the assumption P r = P r (q s , t J , Q r = Q r ( t ) for a par- 
ticular problem. Employing the former assumption rather than the latter (a sort of dual situation) 
is somewhat analogous to the use of the Legendre transformation in partial differential equations 
(Reference 7, V. II, p. 26). 


THE n-BODY PROBLEM 


The rectangular coordinates are z 1? z 2 , z 3 and the corresponding velocity components z 1? z 2 , 
Z 3 ; d\j/dz T is the negative of the force per unit mass due to the self- gravitation of the particles 
themselves. F r (z s , t) is the external force per unit mass (due, for example, to the mass of a 
body not taken as a part of the n-body system and regarded as not affected by the masses of the 
n bodies). Liouville's Equation (Equation 1) takes the form 


3 



r = 1 


3 



0 . 


( 5 ) 


The density distribution function that describes the n-body system is given by 


f 



( 6 ) 
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In Equation 6 z f H and are the position and velocity coordinates ol the a th particle, and z r 
and Z r are the running position and velocity coordinates throughout the space. We shall find the 
form Liouville T s Equation takes under the substitution of Equation 6 with the assumption that z® 
and are functions of all the variables z r , Z r , t, although in practice they would most likely be 
restricted in this regard, often being functions of t alone. A case of some importance (already 
mentioned) is that of taking z as a function of t alone and Z® as a function of both z s and t. 
Examples of this will be given shortly. 

To obtain the density distribution, n, in position space the phase space density function, f , is 
integrated over all velocity space (the concept of density here is the same as that for hydrodynamics): 


N 


00 

JJJf dZj dz 2 dZ 

-00 


3 ' 


Consequently, the potential u due to the n bodies is 


U = - G 


N ( z r , t) dz 1 dz 2 d'Z ' 3 


JJJ /( Z l ~ Z l ) 2 + ( Z 2 “ Z 2) 2 + ( Z 3 " Z 3 ) 2 


co oo n 


= - G 


JJ 

- 00 


£ m„ s (z r - Z r [a] , Z r - z r [a] ) dZj dZ 2 dz 3 dz 1 dx 2 d^ 3 

yu- 


m a S ( 7 r“V) dz l d7 2 dz 3 

a=l 

]/£( Z s -^s ) 2 


/-> ]/£( z s -" z s [a] ) 


In this final expression for u the functions z s ^ no longer involve Z r and z r (if they ever did) since 
these have been removed by integration. 
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Again the following notation is employed: 


<98 1 

( Z r " 


-2H) 

d\ 

('. -*. h ; 

1 


S 


H 


3+s 


c98 1 

( Z r‘ 


-2 r H) 

*1 
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Special attention must be paid to the term (<9U/<9z r ) S s ^ since the delta function changes z s into 
and d\J/dz r is not well defined for z s = z s @ . Actually, the particular terms in (<9u/dz r )8 s ^ 
where this difficulty arises should be taken to be zero. The discontinuity in <9U/<?z r is such that the 
limiting value of this function at the point (z^ , z 2 ®, z 3 ^) in the troublesome term depends on 
the direction of approach. For one set of directions the limiting value fails to be well defined, but 
for the remaining directions it is well defined and turns out to be zero. Indeed, if the factor 
S^z x - is left until last in the consideration of the product, then the result is 
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with the meaning unclear. However, if the product is taken in another order a definite value is 
obtained: 


(,,-*,«)*(, -*,H) 





~N 


3/2 

J 


iU, -tfi) ’.(*,-*») 




13/2 


5 ( Z 1 “Zi [a] )U(z 2 §(23 “ Z 3 [a] ) “ 0 


Moreover, the term physically represents the action of the self- gravitation force of the particle 
and thus plays no role in determining the motion of the system. Consequently, wherever it oc- 
curs it is regarded as zero, and 



/3/a 


- G 



3( 


[a] r \y 
Z LJ “ Z 


E 3 ]^ [ a ] - £ 13) 




In these sums, as indicated, the terms where /3 = a are to be deleted. This is a manifestation of 
the vanishing of the weak functions treated above. The symbol 8 r s in the last term stands for the 
Kronecker delta (not the Dirac delta function) and is unity for r = s and zero for r / s. 
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On the basis of these results Liouville's Equation in rectangular coordinates for an n -body 
system is 




8 



In the event z r @ and Z® depend only on t, Equation 7 becomes 



Inasmuch as the coefficients of 8 s @ and S 3 ^ s are independent of z r and Z r , the left member can 
vanish only if the coefficients of the delta functions are all zero. But this means that 


n 



#CL 


the standard equations for the n-body problem. 

Numerous problems are treated more advantageously in a coordinate system other than the 
rectangular. Several examples will be considered here in which cylindrical coordinates are most 
natural. Let these coordinates be r , 0, z and the components of velocity in these coordinate di- 
rections be r , 0 , z . Liouville's Equation in cylindrical coordinates has the form (Reference 2, 
p. 187) 


di „ di © di „ di /© 2 dV \ di I F0 UU / <?U ^ \ di _ „ /q\ 

ft + 7M + z Fz + - ?7 +F rjaR + ( r r" dV +F a) W + dJ + F z ) Jz ~ 0 • (9) 

The density distribution function representing the n-body system is 

n 

f = J2 m a 8 ( r . z . R-E^, ®-© [ol] , z-z@) . 


The self- gravitation potential is 


u = 


n 

> ]/r 2 + r 


cos + { Z -?®Y 
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Example-Two-Body Problem, Circular Orbits 

In this case n = 2 and the external force is zero. Let 

rM = r x , = r 2 , gfM = 0 ® (t ) , 0 ® = ( t ) , 

gj[i] = ^ ' ©b] = wr , = z-H = rW = ftH = "zU = 2fe) 


0 ( 


and take r 1? r 2 , oj to be constants. This means we are assuming a density distribution function of the form 

2 

f = 2^ m a 5 (r - r a , (9 (t), z, R, © -cor, z) . 

a=l 

Substituting these into Equation 9 and making use of the relation g(x) 8(x - a) - g(a) S(x - a) gives 


66 ® r.i 

- m i s r s 2 [1] - 


m 2 dt 


S 2 ^ + m x a j 5 2 ® + m 2 oj§ 2 ^ + m x ^o; 2 S 


+ m 0 \co z r 


( 2 <™\ s [ 2 ] mi c [l] m2 o [2l ™ , [i] 

V" r 2 - JJ) V J - 77 JO h 5 - 77 57 V J - m i 57 s 6 


* U = 0 


■2 dr) 0 4 J r \ dO °S~ r ^ JB 0 s" J 111 1 dz °6~~ m 2 dz 0 s 

The subscripts on S again stand for partial derivatives. The potential function terms are 

dO n r l - r 3 cos -&&]) 


¥- s [l] 

dr °4 


Gm~ 


t, 


2 + r 2 2 -2r 1 r 2 cos 


§>P - °« 


[ 0 ® - 0 [ 2 ])] 
, cos ( 0 ® - 0 ®) 


5 W 

3/2 °4 


_ g [2] 

1 [r/ + t 2 ~2t 2 t 1 cos (2fH - 2fH)] 3/2 4 


— s M 

d6 S s 




Gm 


r. r 2 sin (&H ~d®) 

- g [1] 

2 r 2 1 2 o 3/2 ® 

+ r 2 - 2r . r 2 cos ^ L J ~ C 


t.* 


r o r 


Gm 


x sin ^ - '8 kJ ) 


[r j 2 + r 2 2 - 2r 2 r x cos - 0®)] 


§ H 

3/2 s 


( 10 ) 


^ s M = |M s t>] 

dz 0 6 dz 6 


Consequently, Equation 10 takes the form 


60® 


" m i ~6T +m 2 " 


+ m 9 r,'Gm 


*)».« , (- 


68® 

»2 ~dF +m 2 aJ 


S + m. s co 2 r , - G 


r, “r« cos 


(flW -fifH) 


1 ^ 2 r 


U H 

3/2 f ° 4 


r 2 - ri cos(fifH - 0 ®) \ 

1 [r7w"-2r 2 r lC os(^-0W)] 3/2 J 


- irij m 2 G 


- m 2 m x G 


, sin («H-«W) 


jr / + r 2 2 - 2r j r 2 cos (& W - $ IsO )J 

{y® -y®) 

(?H -e®)] 3/2 Ss “ 


2 + r 2 2 “ 2 tj r 2 cos 


[r 2 2 + ri 2 -2r 2 r , cos { 0 ® - £W)] 


3/2 S 5 W = 0 


(ID 
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The delta function manipulations used in deriving Equation 11 reduced it to a point where the coef- 
ficients of all the delta functions are independent of r, 0, z, R, 0 , Z; and so these terms are essen- 
tially independent of each other. Consequently, in order that Equation 11 be satisfied these 
coefficients must vanish. From the S and S 2 ® terms it is apparent that 

6 H = cot + Cj , 

6 ^ = cot + c 2 . 

However, from the S s terms we see it is necessary that 0^ - 6^ - 0 or n . Let us choose v 
and take = cot , +77 . The vanishing of the S 4 coefficients after these substitutions gives 


oo* r . - Gm, 


( r i +r 2 Y 


0 , 


co 2 r 0 - Gm 


1 ( r l +r 2)' 


= 0 


From this it is deduced that m 2 r 2 = mi and that 

G m . 

= . rr : 

r 2 ( r l + r 2 ) 


( r l + r 2 ) 2 


The expressions for co can be combined into a symmetric form: 




'G l^i r l +m 2 r ; 


^ r l r 2 ( r 1 + r 2 ) 


Thus we have obtained, by the Liouville Equation approach, the well-known result of the existence 
of a solution to the two-body problem in the form of two particles revolving in circular orbits 
about their common center of gravity. 

Example- Restricted Three-Body Problem 

In this instance n = 3. Take r ^ , 0® , etc. for a = 1 and 2 to be as they were found in the 
previous example. We shall try to get the straight line ("syzygy”) solution. Consequently, let us 
take 7® = r 3 (constant), 0^ = cot, @k] = i to, and z® = r@ = = 0 . Also f = f 0 + f 1 

where f 0 = £ m a S ^ and f t = m 3 5 ^ . This partitions the density distribution function into a 
part, f 0 , governing the motion of the two principal bodies (which already has been found to satisfy 
Liouville’s Equation) and a part, f t , connected with the motion of the third particle. The f t part 
has to be determined, and Liouville's Equation (Equation 9) becomes an equation in i 19 somewhat 
simplified by the fact that f 0 itself satisfies Equation 9. However, despite its appearance Equation 9 
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is nonlinear because of the presence of f in u. This may be indicated by writing u = U(f) = u(f 0 + fi). 
But the functional relation of u on f is essentially that of integration, a linear operation. This 
means u(f 0 + f = u (f 0 ) +u ( f i) . The equation for f x is found to be 


<9fi 

at 


+ R 


dr 


0^ 

+ r dO 


+ Z 


dz 


@2 

r 


w(fj 


dr 


af 1 au(fj) df 


aR 


dr 


o 

aR 


au ( f i) 

dr aR 


R0 

r 


i w(f 0 y 

r W5 


df, 

<50 


1 W(fi) *£o 
r dd a© 


au(f t ) df, au(f,) df 0 


de 


a© 


a z 


az 


au(f 0 ) at, au(fj) at, 

dz dZ dz az 


= 0 


(12) 


The delta function representations for f 0 and f j reduce this to 


r , + r 


| 3 G m l m 3 


13 ) S 4® ” Gm i m 3 


J 3 ^ 4^ Gm 2 m 3 


r 9 + r 


73 § 4 [2] 


0 . 

(13) 


Equating the coefficient of S 4 ® to zero gives the classical equation for r 3 for the straight line 
solution to the restricted three-body problem: 


CO 


2 


r 3 


G m 2 


r 3 + r 2 




0 . 


Note that the last two terms in Equation 13, the S 4 W and S 4 ® terms, cannot be made to vanish. 
This implies that the assumptions on the make-up of f 0 and f were too restrictive to permit an 
exact solution to the three-body problem. The solution must be regarded as approximate, but in 
what sense it is valid is not immediately apparent from Equation 13. A modified problem can be 
formulated for which f x is an exact solution. This is done by assuming the masses of the first 
two bodies so great in comparison with the mass of the third that its effect on the motion of the 
first two is negligible (the "restricted three-body problem"). In the Liouville Equation approach 
this is equivalent to regarding u(f 0 ) as a given external potential and setting the explicit f 0 in 
Equation 12 (i.e. df 0 /<9R, df 0 /d®,df 0 /dZ ) to zero. That is to say, f is merely the density distribu- 
tion function for one particle moving in an external force field set up by the two unaffected 
revolving bodies. The worrisome last two terms in Equation 13 do not appear, whereas the first 
remains unaltered. 


IMPROVEMENT ON APPROXIMATE ORBIT 
Expansions in Terms of Delta Functions 

This section discusses a method for determining corrections to a given approximate orbit of a 
particle or system of particles in an external force field. The procedure is based on the 
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representation of the density distribution function for the particles in the form of a series of delta 
functions and their derivatives. This turns out to be essentially equivalent to employing a power 
series in the differences between the true and approximate values of the coordinates and velocity 
components. 

Suppose the approximate orbit of a single particle is given by z r = z r (t) and z r = z r (t), where 
z r (t) and z r (t) are known. Suppose further that the actual solution is representable in the form 


f 



r 1 r 9 r 3 s 1 S 9 S 


(t)8 l 




3) 


(14) 


where 


g( r l r 2 r 3 S l S 2 S 3) 


= 8 


( r l) 


(*1 - *i) (« a - S 2 ) «<*»> (2, - £3) 8«(z i -Z,) S<**(Z 2 - 2 2 ) 8<W(Z 3 -z 3 ) 


C ( r 3) 


^ 2 ) I 


and the superscripts indicate the v x derivative, etc. (Note the different use of the superscript on 
the delta functions here than in the preceding section.) On occasion it will be convenient to 
use the abbreviation f = J2 c r (t) « Let z r (t) andZ r (t) be coordinates and veloc- 

ity components for the exact solution so that the exact solution to Liouville's Equation 
is f = m 8 (z r - z r , Z r -Z r ) . Thus, we shall be working with the development of this delta func- 
tion for the point (z r , z r ) in the Equation 14 series of derivatives of delta functions about 
the point (z r , Z r ) . 

Expansions in series of the Equation 14 type bear certain formal similarities to expansions 
in terms of orthogonal functions. In particular, the procedure for determining the coefficients in 
the expansion of a given function is essentially the same in each case— multiply by certain func- 
tions and integrate. In expansions like that for Equation 14 the multiplication functions are power 
functions and so the coefficients are determined by "taking moments". 

Suppose we wish to expand a given function g(z r , z r ): 


g( Z r- Zr) 



( Z 1 


= l) 


E ( S 3> 


( Z 3-Z 3 ) 


(15) 


c r s is determined in terms of g by multiplying: 


( Z 1 - Z l) hl ( Z 2 - Z 2) h2 ( Z 3 - Z 3) h3 ( Z l - Z l) kl ( Z 2 " Z 2) k2 ( Z 3 * ( Z r . Z r) 


Z c ^ i ( z 1 - i 1 ) hl ••• (Za-Z^S^ K-0 ••• 5 (S3) ( Z 3- Z 3) 
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Then the following integration is performed: 

00 

JJ-J K --if 1 ■ • • (Z, -Z 3 f 3 Zr) *1 • • ■ ■ ^3 

— OD 


£ % Si { (-1 --if ls(r,) ( z i - *i) • " J (Z 3 -Z3f 3S& ( Z 3 "Zj) dZ 3 . (16) 


The repeated integrals on the right can be evaluated through the formula (Reference 5, p. 10): 

(z - z) h S (r) (z ~ z) = £(z-z) h ] _ S (r > - rjh(z - z) h_1 ] __ S (r_1) + (2 ) [ h (h “ 1 ) ( z ” ^) h ~ 2 ] _ S (r_2) 

- (3) [h(h - 1) (h- 2) (z - z) h “ 3 ] __ S (r_3) + • • • + (-l) h (h)h! S< r “ h > 

- (-l) h (2)h! S< r “ h > . 

This is to be zero if h> r, which is indeed reflected by the fact that (^) , the binomial coefficient, 
should be defined as being zero for h > r . 

Consequently, upon integration 

/>CO /»CO 

I (z - z) h ( z ” z) dz = (“l) h (h)h! Sh" h )(z-z)dz 

•' — CO ^-00 


= c-i) h ( h) h! V ■ 

J co 

S(Z “ z) dz = 1, 

I" S(z - z) dz = f°° S "(z - z ) dz = * • • = 0. Thereupon the right side of Equation 16 becomes 

J — 00 •'-00 


X 1 


Z_, 


• k » c 

3 ' h 1 h 2 h 3 k 1 k 2 k 3 * 


= -■(:;) -C;) h,- 

This supplies the necessary formula for the coefficients in the expansion (Equation 15), 

00 

+l<3 hj ! — k 3 ! JJ-| ( Z 1 - -if' ( Z 3 -Z 3 f 3 g(-r' Z r ) dz J ••• 


(- 1 ) 


dZ, , 


(17) 


and this represents the inversion of Equation 15. 
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This will now be applied to the case at hand; i.e., g(z r , Z r ) = m 8 (z r - z r , Z r -Z r ) and z r = z rJ 
Z r = Z r . The differences between the true and approximate values are 




Z - Z 


(18) 


The coefficients, as given by Equation 17, specialize considerably: 


h. + "*+k, 1 

CD' h, ! • • r k 3 - ! 


( Z 1 " Z 1 +e i) h ‘ ••• ( Z 3~ Z 3 + '7 3 ) k3m6 ( Z r _Z r' Zr‘Z r ) ! •" dZ 3 


hj+— +k 3 poa 

! ! • ■ ■ k 3 T ( z i”^i + t i) 0 ( z i ~ z j j dz j (Z 3 - Z 3 + V 3 ) 3 $ (Z 3 Z 3 ) dZ 3 

J ^ -CO J-CO 


But (z - z + l ) h (z - z) t h >(z - z) , implying that the formula for the coefficients in the expansion 
of m v (z r - z r , z r - Z r ) in terms of the differences e r and 77 r is 


(t) 


m (-1) 

h l ! • 




(19) 


An alternative way of expressing the derivation of these formulas for the coefficients c h k 
is more in the spirit of the Schwartz distribution theory. This method uses the convolution prod- 
uct, f * g (Reference 5, p. 31), instead of the rather loose procedure with integrations involving the 
delta function but devoid of testing functions. In this more satisfactory presentation we have, 
starting with Equation 15: 

( 2 1 “ Z l) h ‘ ( Z 3 " Z 3 ) k3 g( Z r' Z r) 


= 2 Z ( Z > - SO ” 1 K - -I)" - ( Z 3 - Z 3 f 3 (Z, - Z 3 ) 


= E ( " 1)h ‘fc) h > ! S<rrhl) ( Z i- 2 0 ! S (Srt) (Z 3 -2 3 ) 
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Upon taking the convolution with the unit constant function, 

[( z l ~ Z 1 f 1 •" ( Z 3 -Zsf 3 g( Z r- Z r)] 



■ X 


k 3 • X 


'“h . ,k . 


h,+ 

(- 1 ) 1 


hi ! 


k 3 ! 


( 20 ) 


Since the convolution of two functions is well defined if one of them has compact support (such as 
the delta function or any of its derivatives does), the left member of Equation 20 has meaning if 
g(z r , Z r ) has compact support. The general formula from Equation 20 for the coefficients is, of 
course, 


C h.,k l hj ! ■ • • k 3 ! [( z l z l) 1 ( Z 3 ^ 3 ) 3 f ( Z r’ Z r)] * 1 

And if g(z r , z r ) = m 8 (z r -Z r , z r -z r ) then 


[( Z l - Z l) hl •••( Z 3 ~Z 3 ) hmS { Z T ~ Z r- Z r _Z r)] * 1 

(z, - ■■■ (Z 3 -Z 3 ) k3 m [s(z r - z r> Z r - Z r ) * l] 

m ( z l _Z l) hl (^3 ~Z 3 f 3 

hi ko 

- m € x . 

Thus we again arrive at Equation 19. 

These coefficients as given by Equation 19 are of special interest since they are so simply 
related to the errors. Those actually pertinent are the ones for which h x + **• +k 3 = 1, since 
Equation 19 shows that c 100000 = ~m€ lf c 010000 = - me 2 , etc. 


( Z 1 “ z l) hl,, '( Z 3~ Z 3) k3mS ( Z r 


u 


Application to the Equation of Continuity 

To impose the condition of the Equation of Continuity we substitute Equation 14 into Equation 
5 and get, initially, 


ZX - L I] k - 1 ] L k 

r.,s. r.,s. p r.c. p 

i’ i i ’ i ^ i’ i K 

* IX. IX*--**-* 


o . 


(21) 
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(Since we are dealing with the motion of a particle in an external force field, the dU/dz r term in 
Equation 5 does not enter.) The guiding principle behind handling Equation 21 is to reduce the left 
side to a sum of essentially independent terms. This is accomplished by making the coefficients 
of the delta functions independent of the phase space variables z r and Z r . Manipulations that bring 
this about are of the type used earlier. In particular, 




M-.+l.*" 


>(*, -*i) ( Z i-Zx) (Z 2 -Z 2 ) (Z3-Z3) 




s k S 




V 


(ri,~,s k +l,"\s 3 ) 




gMl) 5(^2) g? 3 "^) g( S l , *" ,S k +1, "”> 




i+A 2 +A 3 


dz ^ 2 dz ^ 3 


(evaluated at z r - z ) . 


Employing these relations and, where necessary, shifting the indices on the various sums so that 
all have the delta function derivatives written as enables us to add all these sums as one 

sum on ^ r ‘ ,s ^ with coefficients independent of z r and Z r . These coefficients must vanish if the 
equation is to be satisfied, and so the governing equation for c r s (t) is 


3 




3 

L 

p-i 


( s p + J ) 


-i,***, s+i: 


3 03 on CD 



= 0 . ( 22 ) 

In Equation 22 any c r s with a negative index is regarded as zero. In a particular problem in 
practice it may happen that z p = z p , which simplifies the equation somewhat. The dot over the 
letter, as customary, stands for the total derivative with respect to the time variable. 


An inspection of Equation 22 shows that ultimately an integration will have to be effected in 
determining c r s . In anticipation of this we turn our attention to the constants of integration. 
These depend on the initial values of the differences of z r (t) and z r (t), and Z r (t) and 2 r (t); i.e., 
on e r (0) and i? r (0). Consequently, with these numbers specified the initial values of c r s (t) 
follow immediately from Equation 19: 


1 


S. Si (0) 


_ m(-l) 


„ r i 


(0) £ 2 ^ (0) 6 3 r3 (0) (0) 7)l 2 (0) (0) 


(23) 
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In the special case that is apt to occur, namely that the approximate solution (z r , z r ) has the 
same initial values as the exact solution (z r , z r ) there is considerable simplification: c 000000 (0) = m, 
and all other c. (0) = 0. 

i 1 s i 

Equation 22 is a differential-difference- sum equation in c r s and as such can be studied from 

i * i 

many viewpoints. Perhaps it is most natural to regard it primarily as a difference equation which 
may be treated as a recursion relation, each step in the sequence providing a differential equation. 
A logical successive approximation scheme based on this view begins by considering only the 
terms in Equation 22 which start with c r s whose subscripts add up to unity, r l +r 2 + r 3 + s 1 +s 2 ^ s 3 = 1. 
All terms having indices on the c’s adding up to more than unity are suppressed. As apparent 
from Equation 19, terms of the second and higher degree are neglected in the correction functions; 
i.e., it is a first order theory in the size of the error functions. To obtain a second order theory 
we place the values of c r with + • • * + s 3 = 1 (found by integrating the system of equations 
just indicated) into the equations extracted from Equation 22 headed c r s with r x + • • • + s 3 = 2. 

All terms where this sum is greater than 2 are removed. This system of differential equations is 
solved for c r s (with r 2 + • • * + s 3 = 2) and these values are used in the equations beginning 
c r with r x + • • • + s 3 = 1. Solving this system gives c r s with r x + • • * + s 3 = 1 to the second 

i 1 i i ’ i 

order in the correction functions. This process evidently can be continued in theory with the set 
of equations where the sum is 3, etc. 

Of course, other procedures can be devised for obtaining approximate solutions to Equation 
22. We shall give one that draws heavily on the special form of the equation. The one- dimensional 
problem is treated and then extended to three dimensions. Let 

G r ,A Z ~ F (\ = 0) , 

= (-^ +1 ( r \^) (* = 1 , 2 , 3 , •••) . 

A 

The one-dimensional equation with this simplification in notation is 


- (z-Z)c r _ ls - (s + 1) 


’ r~ 1, s +1 


co 

L 

A=o 


r , A 


r+A , s~l 


(24) 


This can be regarded as a sum equation in c r+ ^ rl , the left side regarded as known by assuming, 
for example, approximate values. In that event this equation is of the "semireduced" or "semi- 
normal M type (Reference 8, p. 535). However, we shall examine Equation 24 from another point 
of view. One of the other ways is to integrate formally and use successive approximations. This 
in effect inverts the d/dt operator. Convergence properties are apt to be improved, however, if 
a larger portion of the differential- difference operator that transforms c r s on the left is inverted. 
This will be done for the operator t where Te r s = c r s - (s + 1) c r _ ls+1 . Consequently, we 
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write 


( S + 1) C r-l,s + l - 'J’r.sOO 


(25) 


where 


•r.s(t) ■ (2-Z> C r-,, S + L • 

A=0 

The simplifying feature of Equation 25 is the fact that the sum of the indices on each of the 
two terms on the left is the same, namely r + s . From the difference equation point of view this 
means that Laplace’s substitution (Reference 8, p. 427) is effective: p - r + s , o- = s. This changes 
Equation 25 to 


C p~cr,ar ( a + ^ ) C J 0-cr-l ,cr+l ^p-CT ,cr ( ^ ) > (26) 

wherein it is apparent that, as a difference equation, this is essentially an equation in a. This 
may be emphasized by writing c p _ c a - u a and = <P a (t). Thus, the equation to be ex- 

amined becomes 


“o- “ (° + !) U a+1 - 0a- (t) . (27) 

The homogeneous equation associated with Equation 27 (i.e., the equation with </> CT (t) replaced 
by zero) has the solution (Reference 8, p. 327) 


T(> + 1) ’ 

where y is an arbitrary constant and j3 is an arbitrary periodic function of a, of period one. Since 
we are concerned with only integral values of a then /3 is a constant. The variation of parameters 
method then suffices in solving inhomogeneous Equation 27. It is assumed that 


U a (t) = r(a+ 1 ) W <7 (t) 


(28) 


and this is substituted into Equation 27. Thus the equation for w a (t) is derived: 


+ r» a 


r(g+i) 

y° e -yt 


<P<T ( * ) • 


(29) 
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The "shifting” operator E, i.e., Ew a = w CT+1 , is introduced and regarded as constant, so that the 
differential- difference equation (Equation 29) takes on the appearance of being only a differential 
equation: 


K ~ y(E- l)w CT = <MO . 


This linear differential equation with constant coefficients has the solution 


-yt 


f erE(t-r) ^ (T) dT + e W-0. C(r 

•/() ' 


(30) 


(31) 


where is independent of t. Although f(o- + 1 ) 7 ^ is independent of r, it is left in the integrand 
since it must lie to the right of e w, ' T ) which operates on it (as shown by the presence of e). Then, 
the solution to Equation 27, by virtue of Equations 28 and 31, has the form 


r 

ro + 1 ) 


f 


e yE(t-r) 1 ) 

7 - 


4> a ( T ) dT + P(cr + 1 ) e7Et C o 


(32) 


The effect of the operator has to be determined to make this usable. This is perhaps 

most easily accomplished by employing the Maclaurin expansion for the exponential function and 
making use of the property E n f(a-) = f(a + n) . We have 


y 0 ' \ ’ 1 , . r(cr + 1 ) 

T(cr + 1) 7 (t- T) k V ^~ 1 'Pa( T ) dT + 

1 k=0 


CO 

k=0 


ryiT) 7 ttt t k E k 


r f “ 

y 17 y -1 1 

r(o- + 1 ) / .in 

J(\ k=0 


v x t T(cr + 1 + k ) 4 y a ^ 1 1 

T ( t - r) k — 0 atk (T) dr + r(cr+ 1) 2_, k 1 


y k t k c a 


I _ 00 03 

( C7 k k )( t " T ) k ^a + k ( T > dT + H> + 1) ^ <’>' t > kC CT +k • 

J Q k=0 k=0 


The symbol y k k j , as before, stands for the binomial coefficient. The quantities C^, which are 
independent of t but dependent on o-, stem from the work on Equation 27. In Equation 26 which in- 
volves both p and <r, p playing no effective role, c +k can be substituted for c a+k . Consequently, 
the general solution to Equation 25, upon retrieving the variables r and s, is found to be under the 
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assumption that term by term integration is justified 


00 00 

: r,s (*) = ( S ^ k ) [ tk * Vk.s-nc (*)] + ry^TT) 


r+s , s+k 


(33) 


For the r integral in this formula the convolution notation has been used for greater succinctness: 


* <E 


r-k , s+k 


(t) 


f 


(t -T) k $ rH( s+k (t) dr 


The initial conditions serve to determine the constants C r s ; for t = 0 in Equation 33 

C r , s (°) = r(s 7 + 1) C r+s , s • 


which implies that 


C 


r + s , s + k 


f( s + k + 1 ) 

yS+k ^r'k,s+k\^J * 


(34) 


These values are incorporated into Equation 33 and the final result is 



If the exact and approximate orbits have the same initial values and the range on r and s is 
taken to be r + s > l , the final series in Equation 35, the one involving c r _ k s+k (0) , drops out com- 
pletely. This is no real restriction since c 00 (t) is the mass and therefore considered known. The 
upper limit on the summation signs in Equation 35 is r instead of co because c r _ k s+k and G r _ k X 
vanish for k > r . 


Although Equation 35 can be said to solve Equation 25 if the right member of Equation 25 is 
taken as given, actually the right side involves the unknowns, c r s . Thus, in effect the inversion 
of the operator has changed the differential-difference equation (Equation 24) into an integral 
equation, but as such it is well adapted for approximation methods, the method of iteration for 
example . 
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For a first order theory the equations are obtained from Equation 35 by setting c r s (t) = 0 
for r + s > 2 . These equations are 


C 01 (t) 


[(^-F) Coo + F z c 10 ] dr + C 01 (0) 


c io (t) 


f (z"Z)c 00 dr+ f (t - r) [(Z-F) Coo + F z cj dr + c 10 (0) + t c 01 (0) 
^0 Jo 


As mentioned above, the terms c 01 (0) , c 10 (0) , and t c 01 (0) disappear if the exact and approximate 
orbits have the same initial values. 


For a second order theory it is assumed c rs (t) = 0 if r + s > 3. The equations in this case 


are 


f [<» 

Jo 


C 01 I K Z c oo + F 2 c 10 F zz c 20 dr + c Q1 (0) 


C 10 (t) 


(z - Z) c oo dr + | (t - t) 


f 


(2 F) c Q0 + F z c 10 + F zz <_ 20 


dr + c io (°) + 1 c oi (°) 


C 02 ( ^ ) 


[(^ -F) c 01 + F z c n ] dr + Cq 2 (0) , 




| [(2-Z)c 01 + (Z-F)c 10 + 2F 2 c 20 ]dr + 2 f (t-r) [(Z - F) c 01 + F z c n ] dr + c n (0) 


+ 2t c 02 (0) _ 


C 20 ( * ) 


f (z-Z)c, 0 d r + f (t-r) [(z-Z)c 01 + (Z-F)c 10 + 2F Z c 20 ~| 
^0 ^0 J 


dr 


J[ (t-^) 2 [(Z"F)c 01 + F,c n ] 


+ (t-^) 2 (Z-F)c 01 + F, Cj j dr + c 20 (0) + tc„ (0) + t 2 c 02 (0) 


The extension from one dimension to three dimensions is almost immediate. After the sub- 
stitutions of r . + s | = p i and s. = o-. a separation of variables is performed as before. In place of 
Equation 25 Equation 22 is now considered, written as 


1*2 '3 3 l a 2 a 3 


(t) 



"l, , " 1 r,,s., "*\s +l p ***,s. 


r l r 2 r 3 s 2 s 2 s 3 


(t) • 
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The result comparable with Equation 35 is 


r, r. 


V V ' \ / s l +k l\/ s 2 +k 2\/ S 3 +k 3\ / k 1 +k 2 +k 3 \ 

r l r 2 r 3 s l s 2 s 3 ^ ^ / . / i / i \ ^1 / \ ^2 /\ ^3 / V ^ r r k l> r 2 -k 2» r 3“ k 3' s l' fk l- s 2 +k 2> s 3 +k 3/ 

kj=0 k,= 0 k 3 =0 


r l r 2 r 3 


E£EE:0CC)CC) 


k i+k 9 +k, 

t 1 " 3 C ( 0 ) 

r r k l' r 2“ k 2‘ r 3“ k 3* s l +k l' s 2 +k 2* s 3 +k 3 


kj=0 k 2 =0 k 3 =0 


The preceding work on the improvement of an approximate orbit for a single particle suffices 
as well for a system of n particles. It is only necessary to let the indices 1, 2, 3 refer to the 
first particle, the indices 4, 5, 6 to the second, and so on, letting the range on p and on i in r i and 
s i be 1 to 3n, instead of 1 to 3. However, it is just as easy to obtain the necessary formulas and 
equations by keeping these ranges 1 to 3 and indicating the different particles by a separate index, 
say a , with range 1 ton, In this latter approach the approximate orbits of the n particles are in- 
dicated by z r = z r ^ (t) , Z r - Z r H (t) where r = 1, 2, 3 and a = l, 2,—, n . The exact solution is 


n co oo oo oo oo oo 

■ EEEZEEE 


r l r 2 r 3 s l s 2 s 3 





b 


( r 3> 


H 





-Z 3 [a] ) '(36) 


for which the abbreviated notation is 


f = 



(t)& 


[a] 


( r 0 



If the actual trajectories are given as z r - z r ^ (t), Z r = (t), then the exact solution also has 

the representation 


n 



The differences between the exact and approximate values for the coordinates and velocity com- 
ponents are e r ^ = z r @ - z r ^ and = Z® -Z r H, and we obtain as a counterpart to Equa- 
tion 19 the formula for the coefficients in Equation 36 for the expansion of the exact solution 
distribution function: 


c 



(t> 


h. + ’^+k, 
(- 1 ) 1 


h, • 


k ’ 
k 3 . 




(37) 
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Working as before leads to the fundamental equation, the counterpart of Equation 22: 



P = 1 A 1 =0 A 2 =0 A 3 =0 Aj A 2 A 3 


= 0 . (38) 

Here as in Equation 22 the notation c ^ >r _ means that the subscripts are the same as for 

c® r c c c ? with the exception of the subscript r o which has been changed to r - 1. Also, as 

r l r 2 r 3 s l s 2 s 3 r-, P p 7 

before, any c r ” a with a negative index is considered zero. 

If the self- gravitation of the particles forms an appreciable part of the force field then this 
self- gravitation effect must be included in the F p ^ term in Equation 38. This self- gravitation 
term has a reasonably simple formulation involving c ^ which will now be derived. 

l’ S i 

The position space density N (z r , t) is obtained by integrating over the velocity space: 


CO 



a r . 


The c r [ a I 2 r 3 ooo C 1 ) involve the masses of the particles, as shown by Equation 37. The gravitation 
potential per unit mass due to the particles themselves is then 


00 



The convention adopted earlier is that the negative of the gradient of the potential is the force 
vector. When -<9U/<9z r is incorporated into F p? (see the earlier treatment of then -body problem) 
for the evaluation of F p at z r = z® , assumed in Equation 38, the terms that do not become well 
defined are suppressed. 
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WEST FORD NEEDLES-DISPENSING 


An example of a problem that can be treated as degenerate in velocity space but not in posi- 
tion space is provided by the West Ford Needles (References 9-13). The intention of the West 
Ford project was to place a large number of small needle-shaped particles into circular orbits 
around the earth and form a thin belt. This presents several interesting dynamics prob- 
lems, in particular that of dispensing the needles from the original container to spread out into 
the belt, and that of the dispersion effect on the belt by the perturbational forces that may be as- 
sumed to be present. It is not the purpose here to attack these problems exhaustively but rather 
to provide a sufficiently broad treatment to exhibit the aptness of the weak functions in situations 
of this general character. 

First the problem of the dispensing of the needles from the container will be treated. The 
natural coordinates are the spherical coordinates p, 9, 0 where 9 is the longitude and 0 the co- 
latitude. The components of the linear velocity vector in these respective directions are P, 0, 
Liouville’s Equation for such a coordinate system is 


d± 

dt 


+ P 


<9 f 


© 


dp p sin 0 


d9 


<D 

+ P 

(- 


if ( 

d4> + \- 


d u 

dp 

i au 

p sinc£ 96 


0 2 
+ 7 

p© 

p 


of. \ii 

P )d P 


cot 0< f>© 

p 


£f ( 

a© + (" 


au 

<90 


Pf 

P 


cot 0© 2 \ <9f 
p ) <9<D 


= 0 


(39) 


We shall deal with an equatorial belt. Thus the problem could be regarded as essentially 
two-dimensional if the belt is assumed to be a curve (as we shall assume). However, the three- 
dimensional equations will be retained and the delta functions will take care of this specialization 
automatically. 

At time t = 0 the needles are expelled from their container at a distance p 0 from the center 
of the earth. They move in the longitudinal direction with a velocity dispersion function, in this 
direction only, of g(@). That is to say, the initial velocity distribution imparted to the particles in 
the 9 direction by the carrier rocket and the expulsion mechanism is denoted by g(@) . This func- 
tion will be taken as normalized so that its integral from -°° to +<» is unity, and it will also be 
taken to be symmetric about the circular orbit velocity p 0 co. Thus, Liouville T s Equation (Equa- 
tion 39) is treated as an initial value problem, the value of f at time t = 0 being 


f(p. 0, <#>, P, ©, <*, 0) = Z(P-P 0 ) 8(6) 5^--5)s(P)S(0)g(®) 


(40) 


In other words the needles, of total mass M, are concentrated initially at the point ( p 0 , 0, 77/2) 
with zero velocity components in the p and0 directions but with a velocity distribution g(©) in the 
£ direction. The Po 2 is present in the denominator because the delta function in curvilinear coordinates, 
when written as the direct product of delta functions of the individual coordinates, must have the 
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Jacobian of the transformation in the denominator (Reference 14, p. 292); for spherical coordinates 
the Jacobian is p 2 sine p which is converted to p Q 2 because of the presence of Z{p- p 0 } and 8(<£-tt/2). 

The method to be used first for this initial value problem is that provided by the Picard proc- 
ess on t, which essentially amounts to obtaining f as a power series in t . This has the advantage 
of being carried out rather easily but it has the disadvantage of presenting f in a form that does 
not display all the properties of the distribution function as clearly as certain other methods, in 
particular the method illustrated later. 

The Picard process consists in writing Equation 39 in the form 


r 

df 0 df <bd± I 

f _5U _0^ d> 2 ' 

\ 11. . / 

1 dU 

p© 

cot ^ 

i di 

l 

dp + p sin 0 66 + p 6<f> + 1 

dp + p + p j 

1 dP + \ 

( p sin 66 

p 

p 1 

1 d® 


+ 


/ ^ 60J 

PO 

cot 00 2 ^ 

i a 

\ P d<t> " 

■ T + 

P ) 

1 d<t> 


dt , (41) 


and employing iteration, starting with f 0 in place of f in the integrand. 

We shall consider the force field due to three sources: the earth's gravitational field, the 
self- gravitational field of the particles themselves, and the gravitational field of the needles' 
original container which presumably will travel with the needles after they have been ejected. 
The potential per unit mass due to the earth's gravitational pull is 


U (1 > 


GM, 

P 


That due to the needles themselves is 



-2-n 

/*00 



N ( p y 6, t ) p 2 s in 4> dp d6 d<£> 

J0 m 

0 «j 

V ' p 2 + p 2 - 2 pp |cos 4>cos + sin 4> sin <j> cos {6 -6) 


where N (p, i 9, <£, t) is the spatial distribution of the needles. On the assumption that the con- 
tainer (with mass M c ) is traveling in a circular orbit in the midst of the needles, its potential 
per unit mass is 


U< 3 > 


GM 

c 

[yO 0 2 + p 2 - 2 p Q p sin 0 cos (0 - cot )] 1/2 
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For the initial approximation f = f 0 , so that 


N(f 0 ) 


f 0 dP d© d4> 


M / TT\ 

— S(P"P 0 ) S (0) S[4>- 2) 8(P) S(<I>) g(0) dP d© d4> 


Pa 


M ( 7T\ 

-j S ( p-po ) 8(0) S (0- “j) 


The initial approximation of the self- gravitation potential consequently is 

/.27T 


U 0 <2) = -G 


IJ 


N Q (p, 0) p 2 sin0dpd0d0 


o ^ o 


)/p 2 + p 2 - 2pp [cos 0 cos 0 + sin0sin0 cos(i9 - <9)J 


GM 


Vp 0 2 + p 2 ~ 2p 0 p sin 0 cos 6 1 

The initial potential due to the container is found by setting t = 0 in the formula for U (3) ; 


u o ( 3 > 


GM 


Vp 0 2 + p 2 - 2p 0 p sin 0 cos B 

Consequently, the combined initial potential is 


G(M + M c ) 


GM* 


yp 0 2 + p 2 - 2p 0 p s in 0 cos 0 

The terms involving the potential, from Equation 41, 


£Up 

ap ap 


c9U 0 af 0 


p sin 0 SB d 0 


JL fUo 
P ^0 (90 


simplify considerably because of the delta functions in f 0 . In fact, 


<?u 0 af 0 


7 1 o GM E M / 77 \ 

ap ap = 77 77 s (^)M 0 -- 2 js'(p)S(o)g( 0 ) 

/^o /°0 ' 


x ^Lo Hs _ i ^Lo TLo 
psin0 a# a® p d 0 aO 


= 0 
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Substituting f 0 for f in the integrand of Equation 41 then leads to the first- order- in - 1 approxima- 
tion to f : 


f (i) = J 7 S(p-p 0 ) 8(fl)s(*--|)8(P)8(*)g(e) 


4 s (r Po) S' (0) s^-f) S(p)5(0)g(©)© 


M ® 2 GMM e\ / 77\ 

+ i— - — ) a(p-Po) 2 ) s ( p ) S(®)g(e) 

S .^0 Po } 


This phase space distribution can be integrated over the velocities to give the position space 
density function through the first power in t. Simplifications occur here because 


I 


S'(P)dP = 0 


and because g(0) is assumed normalized so that its integral is unity. This latter specialization 
also entered into the computation for N 0 . The first order position space density is 

N d> = N ( f O>) = S {P~Po) 2j $'(&) ■ 

The factor p 0 co arises as from the integration of 0g(0), which uses the assumption that g(0) is 
symmetric about the point p 0 co\ 


f 


0g(0) d® 



d® 


" P o w 



g(® ) d® + 



a) g(©) d® 


- Po w ■ 


The last integral vanishes since its integrand is skew- symmetric about p 0 *>; i.e., an odd function 

of ® - p Q co. 

To obtain f (2) this entire procedure is repeated, but f (1) is used in place of f in the integrand 
of Equation 41. This means also, of course, that u (1) = U(f (1) ) is employed. Thus, for the 
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needles 1 self- gravitating potential the potential integral with density N (1) is formed. The evalua- 
tion of this integral depends on a formula used before, a(x) S'(x) = a(0) S'(x) - a'(0) 8(x). The 
first order potential arising from the container's mass is obtained simply by expanding its poten- 
tial function in powers of t and keeping only the first two terms. The complete first order poten- 
tial turns out to be 

GM e G(M + M c ) G(M+M C ) cop Q p sin < p sin 0 

U(1> U ^ f<l) ^ P 1p*+p' - 2p 0 p sin 0 cos ~8 W + P 2 ~ 2 -°o P sin 4> cos e ) 3/2 1 

The details in determining f (2) are similar to those for finding f (1) , although, as expected, 
the work is considerably greater. The resulting formula is 


L (2) 
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Hp-Po) 8(0) s(0-^-)s(P) 8(4.) g(0) 

-JT 8(P-P 0 ) «'(«) 8(*-t) 8 (P) 8 <®)b(®)® 
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The essential problem remains, of course, of extracting the information of interest from this i 

prolix expression. The principal difficulty arises from the meaning of the terms involving the j 

derivatives of the delta functions. The most natural way of investigating these terms is to take i 

the various moments of f . These moments give average values for p 9 p 2 , 0 , 6 2 , etc. for the col- 
lection of particles, and thereby describe the distribution. Different delta derivatives take effect 
for different moments. j 

Of primary interest is the position space density function n ( 2) = N(f (2) ) , generated by inte- J 

grating f (2) over the velocity space. From Equation 42 this is immediately found to be 

i 

N ( 2 ) = 77 S ( p ~ p o) s(p-p 0 ) S'(0) s(*--y)t 

1 f M f C U 2 GM e\ / \ f 7T \ 
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- u 2 g(u) du s(p-p 0 ) S"( 0 ) &(</>- -j) 
*^-00 ' 


2M f U 2 GM E , . / 77 \ 
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(43) 


It was assumed that (@ - p 0 cj) g(0) vanishes at so that the term involving 

g(@) + (®-p 0 w )b'(®) = ik [(®“^o “) g(®)] 


vanishes under integration. 

The integration of n (2) over all position space should, and does, give the total mass of the 
needles: 


fTT r>2 TT /•<» 

I N (2 j p 2 sin 4> dp 6.9 d<£ = M 

Jo Jo Jo 

This provides one check on our work. The masses of the earth and the container do not appear 
since these were not considered to be included in the distribution function. They merely contribute 
the external force field. 

At the initial instant the needles were assumed to be concentrated at one point, but the con- 
tinuous velocity distribution g(@) imposed on them caused them to spread. To get an idea of the 
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extent of this spreading (in the second order theory) in the radial direction, the p - p 0 moment of 
N will be computed. This divided by the mass M gives the average or expected value of p - p 0 at 
time t : 

pTT ~27T *CD 

m (p-Po) 3 (P“Po) n (2) P 2 sin4>dpd& d4> 

•^0 •'0 

= ■ ^(il u2g(u)du -^) 1 p2 ( p ~ p °) S '(p-PoKt* 

1 t 2 , 

^ Pq 

where I is the second moment of g(0) about its point of symmetry, p Q to, 
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since the angular velocity co for a circular orbit is Vgm e p 0 " 3/2 . Therefore the average or expected 
value of p - p 0 at time t is 


(p-Po) = 2^ t2 • 

which, we note, is of the second order in t . 

The spread of the needles in the 9 direction after expulsion from the container is indicated by 
the 6 moments of n ( 2) . It is quickly found that 

_ M P 0 „ 

W - — t - Moot , 

P 0 

so that the average value of 6 is, as expected from the symmetry of the setup, 

6 = cot 
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This, of course, tells nothing about the spread of the needles in the 9 direction. An idea of this 
is gained from the second moment of N (2) with respect to 0; i.e., (when divided by M) the average 
value of 0 2 . The only term that contributes to this in N (2) is the S"(0) term, since 0 2 5(0) = 9 2 s # (0) = 0 
and Q 2 S"(0) = 28(0) . We obtain 


or 


M0 2 


P o 


(l +p 2 a?) t 2 


ife* = i4 2 t 2 + -V it 

' Po 




The moments in the velocity space when computed will give information concerning the aver- 
age values of the different velocity components. 

Our chief dissatisfaction with the power series in t method is that it does not give a compact 
formula for f or N independent of the delta function where such should exist. For example, the 
density distribution in the 6 direction at time t is no longer degenerate but is continuous; i.e. the 
needles have been spread out, and so there should be a function representing this that involves no 
delta functions. Such a function was not obtained by the method just used. A method will now be 
given that does provide this form for the answer. 

Equation 42 indicates that, as far as it goes, f involves <f> and$ only through $(4>-n/2) and 
S(®) (which has the effect of keeping the problem two-dimensional), and P only through S(P) and its 
derivatives. This suggests that f be assumed to be of the form 

00 

f = ^2 E r (p. 0, ©, t) B<0 (P)8(4) 8 (<*>-§■) . (44) 

r - 0 

The advantage of this form is offset somewhat by the fact that the resulting equations are non- 
linear if the self- gravitation of the needles is taken into account. Therefore the treatment here 
will be restricted to the case in which the force field results from the earth’s gravitational field 
alone. Consequently, the equation to consider is a specialization of Equation 39: 

ii nil ® J><3f ( gm e © 2 <d 2 \ df 

dt + ** dp + p sin cj> 36 + p d<t> + l~ p 2 + p + p J dP 

( P® cot0<J>0\df ( P<1) cot 4>® 2 \ df 
+ P J d& + \ p + p ) d<t> ~ 0 ' (45) 


32 


The initial condition is again presented by Equation 40 which implies that for the unknown 
functions E r (p, 6, ©, t) 


E 0 ( P , e, ©, 0) 


P o 


Hp-Po) 8(8) g(®) 


E r (p, e, ®, o) = o , 


(r>0) 


( 46 ) 


The differential-difference equation is obtained for E r (p, Q, ©, t) by substituting Equation 44 
into Equation 45, sliding summation indices so that all the sums contain S (r) (P) S($) S(<£-tt/ 2), 
combining into one sum, and setting the coefficient of S (r) (P) S(4>) S(<£-77/2) to zero. This coef- 
ficient has to vanish in order that the series vanish, since the coefficient is independent of P, <t>, 
<£. This fundamental equation for E r (p, e, ©, t) is 
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(47) 


E r is zero for r < 0. 

This equation can be regarded as a recursion formula in conjunction with an iteration pro- 
cedure. For the first approximation we shall take E x = e 2 = e 3 = • • • = 0; the solution that satis- 
fies Equations 46 is 
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where the series is the Fourier representation of the delta function (Reference 5, p. 33). This 
representation is introduced to facilitate the use of the Laplace Transformation in solving for E/ 1 * 
etc. 


To improve on this approximation to E 0 we determine an approximate Ej , using the value of 
E 0 (1) in the equation for E x (setting E 2 = 0), and then use this value of E t in the equation for E 0 . 
The first three equations taken from Equation 47 are 
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E 0 (1) is the solution to the first of Equations 49 (subject to the side condition imposed by 
Equations 46) with E x set to zero. We thus seek E x (1) as a solution to the second of Equations 49 
with E 2 set to zero, E 0 replaced by E 0 (1) , and the initial values of Ej taken to be zero in agreement 
with Equations 46. That is, the following equation is to be solved for an initial value of zero: 
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The Laplace Transform of Equation 50 with respect to t is 
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By assuming a Fourier expansion for E/ l) , 
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from Equation 51 
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Thus upon taking the inverse transform 
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This same technique is used to determine E 0 < 2 > from the first of Equations 49, 
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This successive approximation scheme for Equations 49 is in essence the same as that sug- 
gested for Equation 22 in the section devoted to the improvement of an approximate orbit of a 
single particle or n particles. c r , s with rj + • • * + s 3 = 1 was of greatest interest for that case. 

i i 

In the present case E 0 is of greatest interest, because it is the only E r that remains upon reduc- 
tion to position space after integrating over the velocities. The expression for N, the position 
space density, is simply 


N(f) 


f dP d© d<t 


E 0 (p, 0, 0, t ) d© 8 [<h - -j) 

^ - 00 


When the value for E 0 given by Equation 52 is employed the spatial density is found to be 

N (2) N ( f ( 2 )) 



The interesting thing about this formula, in comparison with the corresponding Equation 43 
of the previous method, is that here the ", t dependence is expressed through the continuous point 
function g(u) instead of through the delta function and its derivatives. Note that this method has 
transferred the continuous initial velocity distribution in 0 to such a distribution in the position 
space coordinate 

Of course, as in the previous method, it is possible to form the various moments with respect 
to N (2) given in Equation 53 and obtain the same results for p - p 0 , 0 9 (P 9 etc. as in the previous 
method. However, the particular form for Equation 53 enables us to get other information of in- 
terest. For example, the longitudinal variation of the needle density is displayed by the formula 
resulting from integrating N (2) in the p and 4> directions: 



WEST FORD NEEDLES— DISPERSION 


After the dispensing mechanism has spread the needles out in a belt around the earth, the 
question arises of the permanence of the belt under the action of the various disturbing forces. 
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Ideally, the unperturbed distribution has the shape of a circular wire of uniform density in the 
equatorial plane. The particles are to have a nonzero velocity component only in the longitudinal 
direction, of a magnitude which provides a circular orbit under the earth's gravitational field 
alone. The distribution function that describes this is 

M f Tf \ 

f = 8(*) . (54) 

The 2tt in the denominator is necessary so that the integral of f over all phase space is simply M, 
the mass of the needles. 

This f represents an equilibrium distribution under the earth's gravitational force alone and 
so it is a steady state solution to Liouville's Equation. At t = 0 the perturbational forces, repre- 
sented through a potential U, are assumed to take effect, modifying the distribution function by the 
addition of a function F which vanishes at t = 0: 


f = f + F . 


This is equivalent to saying that f is to satisfy an initial value problem under the combined earth's 
gravitational and perturbational force fields with an initial value f as in Equation 54. 

Under these circumstances Liouville’s Equation becomes an equation for F, 
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For the problem to remain simply two-dimensional, the disturbing force in the equatorial 
plane must act only parallel to this plane. A way of expressing this, which will be used here, is 
to have S(<£-7r/2)dU/<?<£ = 0. 

Letting 
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and substituting into Equation 55 gives 
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where U is U evaluated for 0 = tt/2 . This equation is to be solved under the condition that F vanish 
at t = 0. Since we are interested primarily in the character of F in position space it is expedient 
to assume an expansion for F of the form 


F(p, e, p, ©, t) = (p,a,t ) s< r)(p) s c o ( Q-, o . ) . (57 ^ 

r = 0 s - 0 

Thus, in finding the position space density the delta functions in Equation 57 are removed by in- 
tegration, so that N is expressed in terms of u rs (p, #, t). We hope to determine u rs (p, <9, t) in a 
form devoid of delta functions of 0 . 

The equation in u r s , after combining Equations 57 and 56 and making some delta function 
manipulations, is 
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Two types of disturbing forces will be considered: solar radiation pressure, and drag due to 
electromagnetic effects (References 9-13). The solar pressure will be taken as a constant force 
acting in the negative x direction, the earth's shadow effect being neglected. The potential then 
has the form Spcos 6 sin0. The electromagnetic drag can be approximated for particles traversing 
a circular or nearcircular orbit by a conservative force with potential Ad , with A a constant. Con- 
sequently, the disturbing potential is 


U - A# + Spcos# sin <£ . 


(59) 
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The functions u rs are to satisfy the system of equations derived from Equation 58 by using 
Equation 59, 
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under the initial condition of the vanishing of all u rs at the initial instant. This problem submits to 
the technique used earlier on similar problems. The first approximation to u 00 is obtained by 
solving the first of Equations 60 with u 01 and u 10 set to zero. This homogeneous equation with zero 
initial data evidently has as its solution u 0 ( 0 15 = 0. For the first approximation to u 01 and u 10 we put 
u oo = u o y> = 0 and setu 02 = u n = u 20 = 0 in the second and third of Equations 60: 
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This system of linear partial differential equations with constant coefficients (as far as the 
differentiation variables are concerned), subject to the side conditions of zero initial data on t 
and periodicity on 9, can be handled by various techniques. A convenient way is that used previ- 
ously: taking the Laplace transform with respect to t, solving the resulting system of total dif- 
ferential equations in 9, and then taking the inverse transform. The result in this case is 
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Therefore the value of is improved by solving the first of Equations 60 after evaluating 
the right side through the use of Equations 62. The formula obtained for u 00 is 
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For a still better approximation and u^ 1 ) could be found from four, five, and six of 

Equations 60. These values and u 0 ( 0 2) would be used in equations two and three to yield u 0 < 2) and u. 


( 2 ) 


Finally these values would be put in the first equation and u 00 solved. 

An analysis of the results would begin with the determination of the position space density, 
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The variation of the density with 6 is indicated by the integration of N in the p and 4> directions, 
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If the radiation pressure effect alone is considered it is noted from Equation 64 that the density 
runs essentially as W2np 0 (l -k cos 9) 9 meaning that initially the density concentration increases 
on the far side from the sun and decreases on the near side. 
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An insight into the radial direction effects of the disturbing forces may be gained from a study 
of the radial moment. Straightforward computation gives 
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(65) 


The average value of p- p 0 as a function of 0 is found by dividing Equation 65 by Equation 64. 
The electromagnetic drag is seen to contribute a uniform shrinkage to the belt with a change in p 
amounting to - A 3p 0 co 2 (cot) 3 + • * • which, for a moderate value of t, is likely to be very small be- 
cause of the probable small size of A (Reference 13). The solar radiation pressure causes a shift 
of the entire orbit away from the sun together with a much lesser shift in the 90 degree leading 
direction to the direction away from the sun. This may be noted from the negative lead-off coef- 
ficients of cos 6 and sin 6 in Equation 65. The shape of the belt changes too— from a circle to nearly 
a limacon, with the nearest point to the earth almost on the line from the earth to the sun. 


EQUALLY SPACED ECHO SATELLITES 

Now we shall discuss the problem of whether a number of artificial satellites (for example, 
Echo balloons) which are initially equally spaced around a common circular orbit will maintain 
this symmetry to a high degree, or tend to bunch in one region, thereby leaving gaps elsewhere. 

As in the case of the needles' belt geocentric spherical coordinates prove advantageous. How- 
ever, now the plane <f> = n / 2 will not be the earth's equatorial plane but rather the initial orbital 
plane of the satellites. The motion of the satellites will be treated from the viewpoint of perturba- 
tion theory, the dominant force field represented by that for a spherically- symmetric earth. Then 
the unperturbed distribution f satisfies Equation 45, 

n 

f = ji J2 S{p-p 0 ) H8- 0 a) S^-f)s(P) 8(0 -p 0 «) 8(*> , (66) 


where m is the mass of each satellite and & a - 2an/n +a>t . Let F represent that portion of the 
distribution arising from the perturbational forces, so that the total distribution function under 
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all forces is 


f = f + F . 


(67) 


7 p> T e , and designate the components of the disturbing force per unit mass. The fundamental 

equation for F, on the basis of Equations 39, 66, and 67, is 
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One method of continuing is to assume an expansion for F involving p, Q> <j>, P, 0, $> only in the 
delta functions and their derivatives: 
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The differential- difference- sum equation for c r “~. r6 (t), obtained by substituting Equation 69 into 
Equation 68, is quite elaborate for the general perturbing force per unit mass (T p , T d , . A less 
formidable equation results from a procedure similar to that employed in the latter part of the 
last section. This is a compromise. F is expanded in terms of delta functions in P, and $, with 
P 9 6 9 and <f> (and t) left in the to-be- determined coefficients. Then 
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The substitution of Equation 70 into Equation 68 leads to the differential- difference equation, 
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where 

A xy 0 z 0 = 1 if X, y, z = 1, 0, 0 ; A^ 0 = 0 if x, y, z / 1, 0, 0 ; 

A 010 = 1 if x, y, z = 0, 1, 0 ; A 010 = 0 if x, y, z / 0, 1, 0 ; 

A 001 = 1 if X, y, z = 0, 0, 1 ; A 001 = 0 if x, y, z ? 0, 0, 1 . 

The problem is considered as an initial value problem. The satellites traverse the same 
circular orbit with constant angular speed until the "initial instant," t = 0, when the perturbing 
forces take effect. Consequently, Equation 71 is to be solved, subject to the initial condition 

u xyz 0, 0) = 0. 

An inversion of a differential-difference operator portion of Equation 71 will not be attempted, 
although terms could be extracted that bear a resemblance to those in Equation 24, where such an 


42 


inversion was accomplished. Nor shall we be satisfied with merely inverting d/dt. Rather, a tech- 
nique will be used that lies between these extremes. The difference operator aspect will be treated 
as a recursion relation. Inasmuch as it is the 0 variation that interests us most d/dt and d/dO will 
be inverted, with the remaining differential operators taken care of through appropriate delta func- 
tion expansions. 

For the first few values of the integers x, y, and z, Equation 71 is 


d U 000 -Po ^ 11 000 d U 100 1 d u 010 1 ^ u 001 2 _ c ot (p 
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In Equation 73 SSS = T p (p 0 , 5 a , tt/ 2) S(p - p 0 ) s(5-5 a ) S (<£ - 77 / 2 ). The common characteristic 
of each of these groups is the sum of the indices for the two lead-off terms. For Equation 72 it is 
0, for Equations 73 it is 1, for Equations 74 it is 2. Equations 19 and 37 indicate that the sum of 
the indices on u xyz shows the order of its size. The equations within any group are also noteworthy 
in that the only differentiation operator present operating on terms of the same order as the 
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lead-off terms is ■ § -~q. This situation prompts us to treat the equations by inverting 
this operator, recalling that the side conditions are thatu xyz = Ofor t = Oand that u xyz is periodic 
in Q with period 2 n. 

In employing the successive approximation scheme previously suggested, first Equation 72 is 
considered alone, with the initial approximation = 0. It is evident that the solu- 

tion to Equation 72, subject to the stated side conditions, is u 0 < 0 °> = 0. This approximation may be 
improved by revising the approximation to u 100 , u 010 , u 001 . This is done by solving Equations 73 
after setting u x ( y 0 2 } = 0, where x + y + z = 2, and using the best value for u 000 available, = 0. 
That is to say, the following equations are solved: 
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The presence of the delta functions on the right foretells their presence in the solution and so it 
is postulated that 


U x ( yV = £ (1)V -‘ 0( *’ t) 8 (^-^o) S (^'^) ‘ 

a 


Inserting this representation for u x ( ^ into Equations 75 gives the following system of equations 
for (1) v 000 : 
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A convenient technique for handling this system is to employ the Laplace transform with respect to 
t(-) and the finite Fourier transform with respect to d (^) . In anticipation of this the Fourier series 
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for the disturbing force components was inserted on the right of Equations 76. In these series 
the coefficients are constants—the Fourier coefficients of T p (p 0 , < 9 , 77/2), etc. The transformed 
equations form a linear algebraic system, 
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The variable r is the Laplace transform variable counterpart to t, and, similarly, r is the Fourier 
transform variable counterpart to 6. 

Upon taking the inverse transforms the solution to Equations 75 is found to be: 
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J 


where j3 r = K r - i V 2 ^ r and / 3 r = \ r + i . 

Now the original value for u 000 may be improved. The expressions given by Equations 78 are 
used in Equation 72 and this is then solved by the technique just employed. However, in the 
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present case derivatives of the delta functions will appear on the right because of the presence of 
terms like du^j/dp. Consequently, it is necessary to assume 


j (D 
J ooo 



(*• t ) 8 (/> - p 0 ) S (e - d a ) 8 (<£ - \ ) + 
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The value as determined by solving Equation 72 is 
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where 
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te 


1 - (1 + ireot) e~ ir 


and A r is the conjugate of A r with respect to ft ; i.e., the sign of ft is changed wherever it occurs 
with A r . This also holds for B r and C r . Note that the same notation already has been used for j3 r : 
P t = A r - i ftp r and = A r + i ftp r . 

A second order approximation to u 000 may be obtained by the same general procedure, starting 
with Equations 74 under the assumption that u xyz = 0 for x + y + z = 3 and using the approximate 
values already obtained for u xyz for x + y + z = 0 and 1. As may be imagined, the labor is con- 
siderably greater than it was in the previous discussion. 
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In the first order work the formula employed for the distribution function is 
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Furthermore, since this is an n-body problem it is apparent that the distribution function must be 
of the form 
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where p a (t) etc. are the exact solution functions. 

The total mass is found by integrating over all phase space. This would be, from Equation 83, 
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If the alternative form for the distribution function, Equation 82, is used the integral over all phase 
space is 
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Consequently, the last integral must vanish, so that by using Equation 79 
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Upon integrating, 
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In the first order theory this means (refer to Equation 80) that, since the 6 a are inherently 
independent, 
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a fact that inspection of Equations 81 verifies. 

The item of present interest, the behavior of 0 a9 is determined from the 6 moment of f . In 
order that just one 6 a will be obtained, instead of their sum, the 0 integral is taken between the 
limits 0 a - e and d a + e , 
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This moment formed by using the alternative representative of Equation 82 turns out to be 
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But the expression in brackets vanishes (see Equations 84 and 85), and so in the first order 
theory 
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In a similar fashion the velocity component in the 0 direction may be found by taking the 0 
moment, 


®a(0 “ Po" “ (1)v 0°i°0° (^a- 0 • ( 8 7) 

Now these results will be used to determine the growth of asymmetry in the original configu- 
ration of n satellites evenly spaced about a circular orbit. Under the perturbational forces the 
spacing becomes irregular. Various measures of this asymmetry could be devised; the following 
is used here: 


Q(t) " YL “ ( f 'a+l ~‘‘a)] 2 ' (^n+1 

a- 1 

which compares the spacing between successive disturbed satellites with that between successive 
equally spaced satellites. Such a measure is independent of rigid rotation and, therefore, inde- 
pendent of the choice of the satellite corresponding to a = 1. 

Equation 86 may be written 
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with C r given in Equations 81 (see also the remarks following Equation 81). Consequently, 
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tion sign can be taken without this value for the index. Then 
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The geometric series 
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vanishes unless r = - s, in which case it is evidently equal to n. Therefore, 
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The special case (the case of interest) will be considered where the disturbing force com- 
ponents T p (p 0 , <9 a , 77/2) and l e ^ 0 , 6 a , 77/2) are even and odd functions of 8 a , respectively. That 
is to say, since their complex Fourier coefficients are indicated by \ r and j± x , respectively, we 
have \_ r - \ r and = -^ r . Furthermore, since the complex Fourier series 



and 





i r ^ 

e 


must represent real- valued functions, the noted symmetries demand \ r to be pure real and /x r to 
be pure imaginary. For convenience we may write /x r = - io- rJ so that 

= \ ‘ i - /2cr r , 

which is pure real since \ t and cr r are pure real. Simple formulas connect these coefficients for 
positive and negative indices, 
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Similarly, from Equations 81 


C_ r = C r * , and 


C * . 
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The asterisk represents the conjugate with respect to i (the complex conjugate), whereas the bar 
represents the conjugate with respect to ^2. 

From these relations a simple form for W r w_ r and therefore for Q may be obtained, 
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It is a straightforward matter to compute W r W_ r from this formula in terms of t and the 
Fourier coefficients and cr t . The result, when incorporated into Equation 90, gives 


Q(t) = 
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It is not immediately apparent from the form of this expression that it vanishes when t = 0, but 
such is the case. 
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The first three terms in the braces in Equation 92 ultimately provide the growth to Q(t). The 
bulk of the contribution is apt to come from the first few coefficients cr n , which are those for the 
& component of the disturbing force. 

Let us examine the case of the earth and a satellite of the Echo type. The potential due to the 
earth’s gravitational field is 

U = [* "J 2 ( 3 cos 2 <£ e - l) - J 3 ^5 (5 cos 3 4> e - 3 cos 4> e ) + ■ ■ , 


where J 2 and J 3 are numerical constants, the force field is the negative of the gradient of U, and 
0 e is measured from the North Pole. Suppose the satellite’s undisturbed orbital plane intersects 
the equatorial plane along the y axis and the orbital plane lies above the positive x axis. Let the 
declination of the sun relative to the earth’s equator be and the sun’s coordinates in the p, 0 , <f> 
system (i.e., the system based on the undisturbed satellite’s orbit) be K, 0 , Then 

cos 0 e = cos 6 s in <f> s in - "2 + S ^ + cos 4> cos + S : j 

With the force due to the solar radiation pressure falling off inversely as the square of the 
distance, the potential is 


V - k ^o 2 + K 2 - 2pK (cos 6 s in 0 sin <£ 0 + cos cos c£ 0 jj 


■ 1/2 


Upon taking the gradient of the two potentials and evaluating at the point (p 0 , 8 a , 77/2) the 
components of the total disturbing force per unit mass are found to be 
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In these formulas the solar radiation terms were expanded in powers of p/K and terms of a 
higher degree than the first were neglected. The equations can be written in the form of (finite) 
Fourier series in 6 a by employing trigonometric identities to convert powers of the trigonometric 
functions to functions of multiples of the angle. For the p - direction component 
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Since cos e a = (e i6a + e _iSa )/ 2 , the first term is X 0 , is half the coefficient of cos $ a9 k 2 is half 
the coefficient of cos 26 a , k 3 is half the coefficient of cos 3 0 a , and = \. r for r = 1, 2, 3. 

Similarly, 


/ 77 \ 

" 

(p<P 2j 



15 


8P 0 S 


GMJ 3 sin 3 (<£ 0 - 2 + s i) - ’ GMJ 3 sin (^ 0 ~ 2 + S i) + ^2 sin A) 


;in 6„ 


3_ 

15 


2Po 4 


8p 0 5 


77 

\ 3k P 0 

- 2 + Sj 

J + 2K 3 

77 

\” 

” 2 + S i 

j sin 36 


s in 26 


which makes it plain that p 0 = 0, p 1 is half the coefficient of sin d a divided by i, etc., and p, T = -p_ T 
for r = 1, 2, 3. 


54 



As a numerical example we shall take <P 0 = 137.3 degrees, = 0 degrees, p 0 = 1-1/4 earth 
radii, and k/K 2 = 4.75 x 10' 3 cm per second per second. These figures approximate those for 
Echo I (1960 i 1). With these data, from the Fourier expansions for and T 0 the values for the 
coefficients are found to be: 


\ = - 0 . 6 x 10~ 6 p 0 co 2 , 

\ 2 = + 420 x 10" 6 p 0 . 

k 3 - - 0.6x 10~ 6 p 0 o ) 2 , 


cTj = + 2 . 8 X 10“ 6 P 0 a? , 

cr 2 - + 280 x 10“ 6 p Q co 2 , 

cr 3 = -0.4x 10~ 6 p 0 co 2 . 


Of immediate note is the fact that the second coefficient in each set dominates the others. Be- 
cause of this, the formula for the measure of asymmetry, Q, to a sufficient approximation, can be 
written 


Q - 10" 8 (- 28 . 8 r s in V2 r - 40 . 7 cos /2 r + 5 . 09 cos 2 /2 r 

+ 20. 4r 2 + 35.6) , 

where a configuration of thirty (n = 30) satellites is assumed and cot = r . This result shows that 
Q is essentially increasing as t 2 but that there is oscillation about this parabolic curve of period 
V2 77 in r. The amplitude of these oscillations is so small, however, that the curve for Q is mono- 
tone increasing. Therefore, this first order theory for this example indicates that the asymmetry 
in the configuration steadily grows. Naturally, for longer periods of time a second order theory 
would be required to show whether this trend continues or reverses. 


(Manuscript received April 6, 1964) 
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